function [states weights mean_state best_state] = particleFilter(source_mix, target_mix, states, stddev, iters)

best_weight = 0;
best_state = states(1, :);

weights = ones(size(states, 1), 1)/size(states, 1);
for i = 1:iters
    % select base states
    bases = chooseBases(weights, size(states, 1));
    states = states(bases, :);
    
    % sample from process model
    states = states + randn(size(states)).*repmat(stddev, size(states, 1), 1);
    
    % evaluate states with respect to observation
    for j = 1:size(states, 1)
        weights(j) = evaluateState(states(j, :), source_mix, target_mix);
    end
    [temp_best_weight temp_best_index] = max(weights);
    if temp_best_weight > best_weight
        best_weight = temp_best_weight;
        best_state = states(temp_best_index, :);
    end
    
    % normalise importance weights
    weights = weights/sum(weights);
    
    mean_state = sum(states.*repmat(weights, 1, size(states, 2)));
end



function weight = evaluateState(state, source_mix, target_mix)

% translation model so far
source_mix = transformGMM(state, source_mix);
% source_mix.centres = source_mix.centres + repmat(state, source_mix.ncentres, 1);

weight = 1/computeL2Distance(source_mix, target_mix);
% weight = 1/gmm_kl_mc(source_mix, target_mix, gmmsamp(source_mix, 1000));
% weight = 1/sum((source_mix.centres-target_mix.centres).^2, 2);
